/* 
 * File:   objetos.h
 * Author: jonk
 *
 * Created on 22 de noviembre de 2011, 18:28
 */

#include "objetos.h"


//#ifndef INTER_H 
//#define       INTER_H﻿


inline void klingemberg(double dist, double u0, double rcut, double s1, double s2, double *frep, double *epoten)
{
    double dummy;
    double smed=0.5*(s1+s2);
    if (dist<smed+rcut)
    {


        dummy=u0*exp(-100*(dist-smed));
        ///Para las fuerzas
        *frep=-0.1875*dummy;
        ///Para el potencial
        *epoten=0.00375*dummy;
    }
    else
    {
        *frep=0.;
        *epoten=0.0;
    }
}

inline void RShiftedLJ(double dist, double sigma, double epsilon, double rcut, double s1, double s2, double *frep, double *epoten)
{
    double dummy;
    double smed=0.5*(s1+s2);
    double r0,roff,frac2,frac6;
    r0=smed-sigma;
    double rc=smed+2.5*sigma/r0;
    if (dist<rc)
    {

        roff=dist-r0;
        frac2 = (sigma/roff)*(sigma/roff);
        frac6 = frac2*frac2*frac2;

        ///Para las fuerzas
        *frep   = -48.0 * epsilon*frac6*(frac6 - 0.5) / (roff);
        ///Para el potencial
        *epoten= 4.0*epsilon*(frac6*frac6-frac6);
    }
    else
    {
        *frep=0.;
        *epoten=0.0;
    }
}

inline void magnetic(double dist, double dz, double u0, double *fpartxy, double *fpartz, double *epoten)
{
    double Co, Co2;
    double dist3= dist*dist*dist;
    double dist4= dist3*dist;
    double u30dis=3.0*u0/dist4;
    ///Calculamos la fuerza entre iii y jjj
    Co=dz/dist;
        Co2=Co*Co;

        *fpartxy= u30dis*(5.0*Co2-1.0);
        *fpartz= u30dis*(5.0*Co2-3.0);
        *epoten = u0*(1.0-3*Co2)/(dist3);
}

inline void idealWallInteraction(particula part, parametros param, double *fpart, double *epoten)
{
    ///Distancia de la particula a las paredes superior e inferior
    ///Pared superior: zpu=+lz/2+a
    ///Pared inferior: zpd=-lz/2-a
    double dzu, dzd;
    double cpotHS=0.00375;
    double cpotWall=cpotHS;
    double cforceHS=0.1875;
    double cforceWall=cforceHS;
    double dummy;

    dzu=param.lbox[2]-(part.pos[2]+0.5*part.diameter);
    dzd=part.pos[2]-0.5*part.diameter;

    if(dzu<(0.5*part.diameter+param.rcut)) {
        ///Interaccion con pared superior
        dummy=param.u0*exp(-100.*(dzu-0.5*part.diameter));
        ///Para las fuerzas
        *fpart=-cforceWall*dummy;
        ///Para el potencial
        *epoten=cpotWall*dummy;
    }
    if(dzd<(part.diameter+param.rcut)){

        ///La particula interactua con la pared inferior
        dummy=param.u0*exp(-50.*(dzd-0.5*part.diameter));
        ///Para las fuerzas
        *fpart=cforceWall*dummy;
        ///Para el potencial
        *epoten=cpotWall*dummy;
    }
}


inline double interaccion(particula *part, parametros param)
{
    ///Las interacciones consideradas son:
    ///    - Interaccion magnetica entre particulas
    ///    - Interaccion repulsiva particula-particula
    ///    - Interaccion repulsiva particula-pared
        ///Si calcStress = true, se calcula la parte correspondiente al potencial
        ///del tensor de esfuerzos microscopicos y se actualiza en el objeto caja
    ///En el calculo del tensor de esfuerzos solamente se tienen en cuenta
    ///la interaccion magnetica y la repulsion particula-particula, pero no
    ///la interaccion particula/caja. COMPROBAR!!!
    double dist,dist2,ePoten;
    double fx,fy,fz;
    double fcomp, pcomp;
    double dx,dy,dz;
    ///Ponemos las fuerzas de las particulas de la caja a 0
    initForces(part,param);
    ///Inicializamos la energia potencial 
    double fHS, eHS;
    double fMagxy, fMagz, fMagE;
    double massi=1.0;
    double massj=1.0;
    double u0;
    double sigma=0.1;
    fMagxy = fMagz = fMagE=0.0;
    ePoten=0.;
    ///Recorremos todos los pares de particulas
    for(int iii=0;iii<param.npart-1;iii++)
    {

        ///Interaccion de la particula iii con las paredes
        ///Pared superior: zpu=+lz/2+a
        ///Pared inferior: zpd=-lz/2-a
#ifdef CONFINED
        idealWallInteraction(part[iii], param, &fcomp, &pcomp);
        part[iii].force[2]+=fcomp;
#endif
        for(int jjj=iii+1;jjj<param.npart;jjj++)
        {
            u0=param.u0*part[iii].mass*part[jjj].mass;
            ///Distancia con condiciones de contorno periodicas en plano xy pero no en z
            dist= distancia(part[iii],part[jjj],&dx,&dy,&dz,param);
            RShiftedLJ(dist, sigma, param.epsilon, param.rcut, part[iii].diameter,part[jjj].diameter,  &fHS, &eHS);
            //klingemberg(dist, u0, param.rcut, part[iii].diameter,part[jjj].diameter, &fHS, &eHS);
#ifdef MAGNETIC
            magnetic(dist,dz, u0, &fMagxy, &fMagz, &fMagE);
#endif
           // cout<<fMagxy<<"  "<<fMagE<<endl;
            ///Calculamos la fuerza entre iii y jjj
            fx= (fMagxy+fHS)*dx/dist;
            fy= (fMagxy+fHS)*dy/dist;
            fz= (fMagz+fHS)*dz/dist;

            if(fabs(fx)>1000.0 || fabs(fy)>1000.0 || fabs(fz)>1000.0)
            {
                cout<<fx<<"  "<<fy<<"  "<<fz<<endl;
            }

            part[iii].force[0]+=fx;
            part[iii].force[1]+=fy;
            part[iii].force[2]+=fz;

            part[jjj].force[0]-=fx;
            part[jjj].force[1]-=fy;
            part[jjj].force[2]-=fz;


            ///Acualizamos la energia potencial
            ePoten+=eHS+fMagE;


        }

    }
        ///Consideramos la interaccion con las paredes de la ultima particula
        ///no considerada en el bucle anterior
#ifdef CONFINED
         idealWallInteraction(part[param.npart-1], param, &fcomp, &pcomp);

        part[param.npart-1].force[2]+=fcomp;
#endif

        ePoten=ePoten/param.npart;
        return ePoten;
}


//#endif /*INTERACCION_H*/
